use "${clean}newspaper_analysis.dta", clear
keep if north == 1 | south == 1

// recoding Alexandria, VA to the South (technically part of DC until 1844)
replace south = 1 if city == "Alexandria"
replace north = 0 if city == "Alexandria"
replace state = "Virginia" if city == "Alexandria"

// dropping DC as a border area
drop if state == "District Of Columbia"

// standardizing variables
egen county_surban = std(county_urban_percent)
egen county_srugged = std(county_ruggedness)
egen county_spostoffice = std(county_postoffice_area)

// keeping the sample consistent
cap drop sample
gen sample = 1
foreach var of varlist $county_npcontrols {
	replace sample = . if `var' == .
}

// creating some variables
local c = 1860

gen period = year > `c'
replace year = year - `c'

	
// matrices to store the results for graphing
matrix define R = J(5, 3, 0)


* regressions
* ------------

xtset pub_id


// northern baseline re model
qui xtreg singular c.year##i.period##ib1.north ${county_npcontrols} if sample == 1, re cluster(pub_id)
	
qui lincom 1.period // intercept change at cutpoint
matrix R[1,1] = `r(estimate)'
matrix R[1,2] = `r(estimate)' + 1.96*`r(se)'
matrix R[1,3] = `r(estimate)' - 1.96*`r(se)'
	
qui lincom 1.period#c.year // slope change at cutpoint
matrix R[4,1] = `r(estimate)'
matrix R[4,2] = `r(estimate)' + 1.96*`r(se)'
matrix R[4,3] = `r(estimate)' - 1.96*`r(se)'
	
local beta_re = `r(estimate)'
local beta_re : display %4.3f `beta_re'
local se_re = `r(se)'
local se_re : display %4.3f `se_re'	
	

// northern baseline fe model
qui xtreg singular c.year##i.period##ib1.north ${county_npcontrols} if sample == 1, fe cluster(pub_id)

qui lincom 1.period // intercept change at cutpoint
matrix R[2,1] = `r(estimate)'
matrix R[2,2] = `r(estimate)' + 1.96*`r(se)'
matrix R[2,3] = `r(estimate)' - 1.96*`r(se)'
	
qui lincom 1.period#c.year // slope change at cutpoint
matrix R[5,1] = `r(estimate)'
matrix R[5,2] = `r(estimate)' + 1.96*`r(se)'
matrix R[5,3] = `r(estimate)' - 1.96*`r(se)'

local beta_fe = `r(estimate)'
local beta_fe : display %4.3f `beta_fe'
local se_fe = `r(se)'
local se_fe : display %4.3f `se_fe'


* graphing 
* --------
svmat R

// temp variable to index estimates
cap drop temp*
gen tempn = _n if R1 != .

gen tempslope = tempn > 3

gen tempfe = mod(tempn, 3)
recode tempfe (1=0) (2=1) (0=.)

 
twoway ///
	(scatter R1 tempn if tempslope == 0 & tempfe != .,	msize(large) mlwidth(none) mcolor("${blue}") 	msym(D) yaxis(1 2) ) ///
 	(rcap R2 R3 tempn if tempslope == 0 & tempfe == 0,   lwidth(thick) lcolor("${blue}") 	lpattern(dash) yaxis(1 2)) ///
 	(rcap R2 R3 tempn if tempslope == 0 & tempfe == 1,   lwidth(thick) lcolor("${blue}") 	lpattern(solid) yaxis(1 2)) ///
	, ///
	yline(0, lcolor(gs2) lpattern(dash) lwidth(thick)) ///
	legend(off) ///
	xlab(1 "Baseline" 2 "Fixed Effects") ///
	xtitle("") 	///
	ylab(none, axis(2)) ///
	ylab(-.1(0.05).1, axis(1) grid glcolor(gs7%50)) ///
	ytitle("") ///
	xscale(range(0.5 2.5)) ///
	subtitle("Intercept Shift: 1860 - 1861", box bexpand bcolor(gs2) color(white) size(medsmall)) ///
	fxsize(100) fysize(100) ///
	name(intercepts, replace) nodraw
	
	
twoway ///
	(scatter R1 tempn if tempslope == 1 & tempfe != .,	msize(large) mlwidth(none) mcolor("${blue}")	msym(S) yaxis(1 2)) ///
	(rcap R2 R3 tempn if tempslope == 1 & tempfe == 0,   lwidth(thick) lcolor("${blue}") 	lpattern(dash) yaxis(1 2)) ///
	(rcap R2 R3 tempn if tempslope == 1 & tempfe == 1,   lwidth(thick) lcolor("${blue}") 	lpattern(solid) yaxis(1 2)) ///
	, ///
	text(0.003 4 "{&beta} = `beta_re'") ///
	text(0.0015 4 "se = `se_re'") ///
	text(0.003 5 "{&beta} = `beta_fe'") ///
	text(0.0015 5 "se = `se_fe'") ///
	yline(0, lcolor(gs2) lpattern(dash) lwidth(thick)) ///
	legend(off) ///
	xlab(4 "Baseline" 5 "Fixed Effects") ///
	xtitle("") 	///
	ylab(none, axis(1)) ///
	ylab(-0.012(0.006)0.012, axis(2) grid glcolor(gs7%50)) ///
	ytitle("") ///
	xscale(range(3.5 5.5)) ///
	subtitle("Slope Change after 1860", box bexpand bcolor(gs2) color(white) size(medsmall)) ///
	fxsize(100) fysize(100) ///
	name(slopes, replace)  nodraw
	
graph combine intercepts slopes,  scale(1.5) xsize(100) ysize(50)
graph export "${output}Fig3_np_north_1860.pdf", as(pdf) replace


